Multi-omics integration with Giotto Suite
Joselyn Chávez
2023-08-02
1 Data
To run this vignette, we will use the Skin dataset from the article High-plex protein and whole transcriptome co-mapping at cellular resolution with spatial CITE-seq
2 Installation
To run this vignette, you must download the last version of Giotto Suite
remotes::install_github("drieslab/Giotto@suite_dev")3 Individual modalities preparation
3.1 Create Giotto object
Load Giotto
library(Giotto)## Giotto Suite 3.3.1instrs <- createGiottoInstructions(save_plot = FALSE,
show_plot = FALSE)Create Giotto object using RNA and Protein expression, as well as spatial positions.
my_giotto_object <- createGiottoObject(
expression = list(rna = list(raw = "data/expression_rna.csv"),
protein = list(raw = "data/expression_protein.csv")),
expression_feat = list("rna", "protein"),
spatial_locs = "data/positions.csv",
instructions = instrs)3.2 Add tissue image
my_giotto_image <- createGiottoImage(gobject = my_giotto_object,
do_manual_adj = TRUE,
scale_factor = 0.0625,
mg_object = "data/skin.jpg",
negative_y = TRUE)## do_manual_adj == TRUE
## Boundaries will be adjusted by given values.my_giotto_object <- addGiottoImage(gobject = my_giotto_object,
images = list(my_giotto_image),
spat_loc_name = "raw")Visualize image
spatPlot2D(my_giotto_object,
show_image = TRUE,
point_size = 2)
3.3 RNA modality
3.3.1 Filtering
my_giotto_object <- filterGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_threshold = 1,
feat_det_in_min_cells = 1,
min_det_feats_per_cell = 1)## completed 1: preparation
## completed 2: subset expression data
## completed 3: subset spatial locations
## completed 4: subset cell (spatial units) and feature IDs
## completed 5: subset cell metadata
## completed 6: subset feature metadata
## completed 7: subset spatial network(s)
## completed 8: subsetted dimension reductions
## completed 9: subsetted nearest network(s)
## completed 10: subsetted spatial enrichment results
##
## Feature type: rna
## Number of cells removed: 0 out of 1691
## Number of feats removed: 363 out of 154863.3.2 Normalization
my_giotto_object <- normalizeGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
norm_methods = "standard",
scalefactor = 10000,
verbose = TRUE)##
## first scale feats and then cells3.3.3 Gene statistics
my_giotto_object <- addStatistics(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna")3.3.4 HVF identification
my_giotto_object <- calculateHVF(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_values = "normalized",
method = "cov_groups",
nr_expression_groups = 20,
zscore_threshold = 1.5)## return_plot = TRUE and return_gobject = TRUE
##
## plot will not be returned to object, but can still be saved with save_plot = TRUE or manually3.3.5 Principal component analysis (PCA)
my_giotto_object <- runPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_values = "normalized",
reduction = "cells",
feats_to_use = "hvf",
name = "rna.pca")## "hvf" was found in the feats metadata information and will be used to select
## highly variable features## class of selected matrix: dgCMatrix
## [1] "finished runPCA_factominer, method == factominer"3.3.6 Uniform manifold approximation projection (UMAP)
my_giotto_object <- runUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
expression_values = "normalized",
reduction = "cells",
dimensions_to_use = 1:10,
dim_reduction_name = "rna.pca")3.3.7 Create nearest network
my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
type = "kNN",
dim_reduction_name = "rna.pca",
name = "rna_kNN.pca",
dimensions_to_use = 1:10,
k = 20)3.3.8 Find Leiden clusters
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
nn_network_to_use = "kNN",
network_name = "rna_kNN.pca",
name = "leiden_clus")3.3.9 Plot PCA
plotPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
dim_reduction_name = "rna.pca",
cell_color = 'leiden_clus',
title = "RNA PCA")
3.3.10 Plot UMAP
plotUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
cell_color = 'leiden_clus',
point_size = 2,
title = "RNA Uniform Manifold Approximation & Projection (UMAP)",
axis_title = 12,
axis_text = 10)
3.3.11 Plot spatial locations by cluster
spatPlot2D(my_giotto_object,
show_image = TRUE,
point_size = 1.5,
cell_color = "leiden_clus",
title = "RNA Leiden clustering")
3.4 Protein modality
3.4.1 Filtering
my_giotto_object <- filterGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_threshold = 1,
feat_det_in_min_cells = 1,
min_det_feats_per_cell = 1)## completed 1: preparation
## completed 2: subset expression data
## completed 3: subset spatial locations
## completed 4: subset cell (spatial units) and feature IDs
## completed 5: subset cell metadata
## completed 6: subset feature metadata
## completed 7: subset spatial network(s)
## completed 8: subsetted dimension reductions
## completed 9: subsetted nearest network(s)
## completed 10: subsetted spatial enrichment results
##
## Feature type: protein
## Number of cells removed: 0 out of 1691
## Number of feats removed: 0 out of 2833.4.2 Normalization
my_giotto_object <- normalizeGiotto(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
scalefactor = 10000,
verbose = T)##
## first scale feats and then cells3.4.3 Gene statistics
my_giotto_object <- addStatistics(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized")3.4.4 Principal component analysis (PCA)
my_giotto_object <- runPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized",
scale_unit = T,
center = F,
method = "factominer")## "hvf" was not found in the gene metadata information, all genes will be used## class of selected matrix: dgCMatrix
## [1] "finished runPCA_factominer, method == factominer"3.4.5 Uniform manifold approximation projection (UMAP)
my_giotto_object <- runUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
expression_values = "normalized",
dimensions_to_use = 1:10)3.4.6 Create nearest network
my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
type = "kNN",
name = "protein_kNN.pca",
dimensions_to_use = 1:10,
k = 20)3.4.7 Find Leiden clusters
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
nn_network_to_use = "kNN",
network_name = "protein_kNN.pca",
name = "leiden_clus")3.4.8 Plot PCA
plotPCA(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
dim_reduction_name = "protein.pca",
cell_color = 'leiden_clus',
title = "Protein PCA")
3.4.9 Plot UMAP
plotUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
cell_color = 'leiden_clus',
dim_reduction_name = "protein.umap",
point_size = 2,
title = "Protein Uniform Manifold Approximation & Projection (UMAP)",
axis_title = 12,
axis_text = 10)
3.4.10 Plot spatial locations by cluster
spatPlot2D(my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
cell_color = "leiden_clus",
point_size = 1.5,
show_image = TRUE,
title = "Protein Leiden clustering")
3.5 Save Giotto object
saveGiotto(my_giotto_object, "multiomics_Giotto_object")4 Multi-omics integration
4.1 Load Giotto object
Download the file from https://drive.google.com/file/d/1GyhKWl14MxEJCw_2-e8Bge3d3zY7XLkc/view?usp=share_link
my_giotto_object <- loadGiotto("multiomics_Giotto_object")4.2 Calculate WNN
my_giotto_object <- runWNN(my_giotto_object,
modality_1 = "rna",
modality_2 = "protein",
pca_name_modality_1 = "rna.pca",
pca_name_modality_2 = "protein.pca",
k = 20,
verbose = TRUE)## The NN network name was not specified, default to the first: "rna_kNN.pca"## The NN network name was not specified, default to the first: "protein_kNN.pca"## [1] "Calculating rna - rna distance"
## [1] "Calculating protein - protein distance"
## [1] "Calculating low dimensional cell-cell distances for rna"
## [1] "Calculating low dimensional cell-cell distances for protein"
## [1] "Calculating within-modality prediction"
## [1] "Calculating cross-modality prediction"
## [1] "Calculating Jaccard similarities"
## [1] "Calculating kernel bandwidths"
## [1] "Calculating modality weights"
## [1] "Calculating WNN"
## [1] "Calculating WNN normalization"
## [1] "Calculating WNN graph"
## [1] "Saving WNN results"4.3 Create integrated UMAP
my_giotto_object <- runIntegratedUMAP(my_giotto_object,
modality1 = "rna",
modality2 = "protein",
spread = 5,
min_dist = 2)## [1] "Calculating integrated Nearest Neighbors"
## [1] "Creating igraph"
## [1] "Calculating integrated UMAP"4.4 Calculate Leiden clusters
my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
nn_network_to_use = "kNN",
network_name = "integrated_kNN",
name = "integrated_leiden_clus",
resolution = 0.4)4.5 Plot integrated UMAP
plotUMAP(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
cell_color = 'integrated_leiden_clus',
dim_reduction_name = "integrated.umap",
point_size = 1.5,
title = "Integrated UMAP using Integrated Leiden clusters",
axis_title = 12,
axis_text = 10)
4.6 Plot integrated spatial locations by cluster
spatPlot2D(my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
cell_color = "integrated_leiden_clus",
point_size = 1.5,
show_image = TRUE,
title = "Integrated Leiden clustering")
5 Interactive selection
5.1 Add tissue image
my_giotto_image <- createGiottoImage(gobject = my_giotto_object,
do_manual_adj = TRUE,
scale_factor = 0.0625,
mg_object = "data/skin.jpg",
negative_y = TRUE)
my_giotto_object <- addGiottoImage(gobject = my_giotto_object,
images = list(my_giotto_image),
spat_loc_name = "raw")The interactive selection tool will allow you to run a Shiny app and plot polygons over the tissue.
Note that you can use the blue bars to Zoom in/out at areas of interest over the tissue
You can choose multiple polygon names, one for each drawn polygon.
When you have finished plotting polygons, click on the Done button
my_polygon_coordinates will store a data.table object with the polygon names and xy coordinates for downstream analysis.
5.2 Create a spatial plot
my_spatPlot <- spatPlot2D(gobject = my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
show_image = TRUE,
cell_color = 'leiden_clus',
point_size = 2,
point_alpha = 0.75)
my_spatPlot
5.3 Run interactive tool
my_polygon_coordinates <- plotInteractivePolygons(my_spatPlot, color = "black")
5.4 Add polygons information to Giotto object
my_giotto_polygons <- createGiottoPolygonsFromDfr(my_polygon_coordinates,
name = 'selections')## Selecting col "name" as poly_ID column## Selecting cols "x" and "y" as x and y respectivelymy_giotto_object <- addGiottoPolygons(gobject = my_giotto_object,
gpolygons = list(my_giotto_polygons))addPolygonCells will modify the cell_metadata slot to indicate if a certain cell is located within a polygon area (polygon1,2,3..) or not (no_polygon)
my_giotto_object <- addPolygonCells(my_giotto_object,
polygon_name = 'selections')##
## These column names were already used: nr_feats perc_feats total_expr
## leiden_clus integrated_leiden_clus
## and will be overwrittenVerify that cell metadata has been modified
head(pDataDT(my_giotto_object))## cell_ID nr_feats perc_feats total_expr leiden_clus integrated_leiden_clus
## 1: 10x14 70 0.4628711 474.8438 2 7
## 2: 10x15 237 1.5671494 1228.7768 4 7
## 3: 10x16 92 0.6083449 599.3736 2 5
## 4: 10x17 87 0.5752827 581.0906 5 4
## 5: 10x18 132 0.8728427 807.0911 7 1
## 6: 10x19 336 2.2217814 1601.3652 3 1
## selections
## 1: polygon 2
## 2: polygon 2
## 3: polygon 2
## 4: polygon 2
## 5: no_polygon
## 6: no_polygon5.5 Retrieve cells from the polygon area
my_polygons_cells is a terra::spatVector object that contains the information of each cell within the polygon areas
my_polygons_cells <- getCellsFromPolygon(my_giotto_object,
polygon_name = 'selections')
my_polygons_cells## class : SpatVector
## geometry : points
## dimensions : 465, 4 (geometries, attributes)
## extent : 1, 50, -43, -1 (xmin, xmax, ymin, ymax)
## coord. ref. :
## names : sdimx sdimy cell_ID poly_ID
## type : <int> <int> <chr> <chr>
## values : 10 -14 10x14 polygon 2
## 10 -15 10x15 polygon 2
## 10 -16 10x16 polygon 25.6 Differential Expression between selected regions
Find marker features for each polygon
scran_results <- findMarkers_one_vs_all(my_giotto_object,
spat_unit = "cell",
feat_type = "rna",
method = "scran",
expression_values = "normalized",
cluster_column = "selections",
min_feats = 10)## using 'Scran' to detect marker feats. If used in published research, please cite:
## Lun ATL, McCarthy DJ, Marioni JC (2016).
## 'A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor.'
## F1000Res., 5, 2122. doi: 10.12688/f1000research.9501.2.##
## start with cluster no_polygon
##
## start with cluster polygon 1
##
## start with cluster polygon 2top_genes <- scran_results[, head(.SD, 2), by = 'cluster']$featscomparePolygonExpression will create a heatmap to observe the differential expression between polygon areas
comparePolygonExpression(my_giotto_object,
selected_feats = "top_genes")## using 'Scran' to detect marker feats. If used in published research, please cite:
## Lun ATL, McCarthy DJ, Marioni JC (2016).
## 'A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor.'
## F1000Res., 5, 2122. doi: 10.12688/f1000research.9501.2.##
## start with cluster no_polygon
##
## start with cluster polygon 1
##
## start with cluster polygon 2
You can isolate each region to plot cell metadata such as the number of features per cell.
spatPlot2D(my_giotto_object,
cell_color = 'nr_feats',
group_by = 'selections',
color_as_factor = FALSE,
point_size = 2,
coord_fix_ratio = 1,
cow_n_col = 2,
show_legend = TRUE)
Or plot gene expression per cell by selected region (polygon).
spatFeatPlot2D(my_giotto_object,
expression_values = 'scaled',
feats = c("APOC1", "TMEM132D"),
group_by = 'selections',
point_size = 2)
5.7 Cell types enriched between selected regions
You can compare the percent of cell types by selected region (polygon).
compareCellAbundance(my_giotto_object,
spat_unit = "cell",
feat_type = "rna")
You can isolate each region and plot cell types
## cell_ID selections
## 1: 10x14 polygon 2
## 2: 10x15 polygon 2
## 3: 10x16 polygon 2
## 4: 10x17 polygon 2
## 5: 10x18 no_polygon
## ---
## 1687: 9x46 no_polygon
## 1688: 9x47 no_polygon
## 1689: 9x48 no_polygon
## 1690: 9x49 no_polygon
## 1691: 9x50 no_polygonspatPlot2D(my_giotto_object,
spat_unit = "cell",
feat_type = "protein",
cell_color = 'integrated_leiden_clus',
group_by = 'selections',
color_as_factor = TRUE,
point_size = 1,
coord_fix_ratio = 1,
cow_n_col = 2,
show_legend = TRUE)
5.8 Plot polygons again
Finally, you can plot the regions previously selected.
plotPolygons(my_giotto_object,
x = my_spatPlot,
polygon_name = 'selections')